function [lb,ub] = fsel_linear(X,Y,A,set0,kmax)
% Forward selection for inference for linear function of coefficients
% Last edit: 3/28/17 CBH
% [lb,ub] = fsel_linear(X,Y,A,set0,K)
% X - design matrix
% Y - outcome
% A - linear combination of form A'beta
% set0 - set of initial regressors
% kmax - number of forward selection steps to take

p = size(X,2);
varind = (1:p)';

supp_lb = set0;
supp_ub = set0;

lb = zeros(kmax,1);
ub = zeros(kmax,1);

for jj = 1:kmax
    varind_lb = setdiff(varind,supp_lb);
    varind_ub = setdiff(varind,supp_ub);
    
    ci_lb = zeros(numel(varind_lb),1);
    ci_ub = zeros(numel(varind_ub),1);
    
    for kk = 1:numel(varind_lb)

        use_lb = [supp_lb;varind_lb(kk)];
        zlb = X(:,use_lb);
        b_lb = zlb\Y;
        e_lb = Y - zlb*b_lb;
        [~,V_lb] = hetero_se(zlb,e_lb,inv(zlb'*zlb));
        BLB = zeros(p,1);
        VLB = zeros(p,p);
        BLB(use_lb) = b_lb;
        VLB(use_lb,use_lb) = V_lb;

        use_ub = [supp_ub;varind_ub(kk)];
        zub = X(:,use_ub);
        b_ub = zub\Y;
        e_ub = Y - zub*b_ub;
        [~,V_ub] = hetero_se(zub,e_ub,inv(zub'*zub));
        BUB = zeros(p,1);
        VUB = zeros(p,p);
        BUB(use_ub) = b_ub;
        VUB(use_ub,use_ub) = V_ub;
        
        ci_lb(kk,1) = A'*BLB - 1.96*sqrt(A'*VLB*A);
        ci_ub(kk,1) = A'*BUB + 1.96*sqrt(A'*VUB*A);
        

    end
    imin = find(ci_lb == min(ci_lb));
    imax = find(ci_ub == max(ci_ub));
    imin = imin(1);
    imax = imax(1);
    
    supp_lb = [supp_lb ; varind_lb(imin)]; %#ok<*AGROW>
    supp_ub = [supp_ub ; varind_ub(imax)];
        
    lb(jj,1) = min(ci_lb);
    ub(jj,1) = max(ci_ub);
end
